R パッケージ一覧library(broom)
library(car)
library(jtools)
library(normtest)
library(patchwork)
library(stargazer)
library(tidyverse)
library(zoo)BLUE: Best Linear Unbiased Estimator) とは、次の 5
つの仮定が満たされるときの最小二乗法による回帰係数の推定量のことBLUE
に即して、社会科学において要求される最小二乗法による回帰分析の前提と回帰モデルの妥当性の診断方法を解説する| 回帰分析の前提 | 妥当性の診断方法 |
|---|---|
| 0. 回帰モデルは妥当なモデルである | 学問上の知識に基づいて判断 |
| 1. 誤差項の期待値がゼロ | |
| 2. 加法性と線形性 | Scatter plotで確認 |
| 3. 誤差が独立である | 慎重に交絡変数を検討する |
| 4. 誤差の分散が均一である | 残差プロットによる診断 |
| 5. 誤差が正規分布に従う | 正規QQプロットによる診断 |
| 6. 完全な多重共線性がない | 分散拡大係数 (VIF) による診断 |
・因果推論するためには仮定 1 と仮定 2
を満たすことが特に重要
・この仮定を完全に満たすことは困難 → 傾向スコアを使った回帰分析
RDD: Regression Discontinuity Design)DID: Difference-in-differences)IV: Instrumental Variables)\[Y_i = \alpha_0 + β_1X_1i +u_i\]
\(u_i\) はこのモデルにおける誤差項
「誤差項の期待値がゼロでない」場合でも、誤差項の期待値がゼロになることを示す
つまり、誤差項の期待値 \(E(u_i) = σ(シグマ)、σ≠0\)
ということ
ここでは、誤差項 \(u_i\) の期待値がゼロでない場合から始めても、誤差項 \(ε_i\)の期待値はゼロとなり仮定は満たされることを示してみる
| モデル | 解説 |
|---|---|
| \(Y_i = \alpha_0 + β_1X_1i +u_i\) | 誤差項 \(u_i\) の期待値がゼロでないと仮定 |
| \(Y_i = \alpha_0 + β_1X_1i +u_i + \sigma - \sigma\) | \(\sigma=0\)なので、計算上右辺に \(\sigma\) を足して引く |
| \(Y_i = (\alpha_0 + \sigma) + β_1X_1i + (u_i - \sigma)\) | 新たな \(Y\) 切片は \((\alpha_0 + \sigma)\)、新たな誤差項は \((u_i - \sigma) になる\) |
| \(\beta_0 = (\alpha_0 + \sigma), ε_i = (u_i - \sigma)\)と定義する | |
| \(Y_i = \beta_0 + β_1X_1i + ε_i\) | 誤差項 \(u_i\) の期待値がゼロでない場合、\(Y\) 切片の値は \(\alpha_0\) から \(\beta_0 = \alpha_0 + σ\) に変わるが、傾きは \(\beta_1\) のまま変化せず |
| \(E[ε_i]=E[u_i-\sigma]\) | 再定義した誤差項 \(ε_i = u_i - \sigma\) の期待値をとる |
| \(E[ε_i]=E[u_i] - E[\sigma]\) | 期待値の加法性 \(E[X+Y] = E[X] + E[Y]\) より |
| \(E[ε_i]=\sigma - \sigma\) | \(E[u_i] = \sigma\)であり、\(E[定数] = 定数\)だから(\(\sigma\) は定数) |
| \(E[ε_i] = 0\) | 従って、誤差項 \(ε_i\)の期待値はゼロとなり、仮定は満たされた |
ポイント
・誤差項 \(u_i\)
の期待値がゼロでない場合、\(Y\) 切片の値は \(\alpha_0\) から \(\beta_0 = \alpha_0 + σ\)
に変わるが、傾きは \(\beta_1\)
のまま変化しない
・新たに定義し直した誤差項 \(ε_i = u_i - \sigma\)
の期待値はゼロと見なすことができる
set.seed(20220225) # seed の設定(数字は任意)
n <- 1000 # 生成するデータ数を n = 1000 と設定
a0 <- 1 # モデルの切片 (α0)を 1 に設定
b1 <- 2 # モデルの傾き (β1) を 2 に設定
x1 <- rnorm(n) # x1 の値をランダムに 1000個生成u1 <- rnorm(n, 0, 1) # 誤差項の期待値を 0 と設定
# 標準偏差を 1 と設定hist(u1)
- \(u1\) の記述統計を確認
summary(u1) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.93364 -0.62634 0.02522 0.03083 0.68300 3.02237
当初指定したとおり、誤差項 \(u1\) の期待値 \(E(u1)\) はほぼゼロ: 0.03083
応答変数 \(y1\) に説明変数 \(x1\) を回帰させて結果を表示させる
y1 <- a0 + b1*x1 + u1 # 応答変数 (y1) を定義する
model_1 <- lm(y1 ~ x1)u2 <- rnorm(n, 10, 1) # 誤差項の期待値を 10 と設定
# 標準偏差を 1 と設定hist(u2)summary(u2) Min. 1st Qu. Median Mean 3rd Qu. Max.
6.688 9.267 9.986 9.976 10.644 12.987
当初指定したとおり、誤差項 \(u2\) の期待値 \(E(u2)\) はほぼ 10: 9.976
応答変数 \(y2\) に説明変数 \(x1\) を回帰させて結果を表示させる
y2 <- a0 + b1*x1 + u2 # 応答変数 (y2) を定義する
model_2 <- lm(y2 ~ x1)stargazer(model_1, model_2,
type = "html",
digits = 2)| Dependent variable: | ||
| y1 | y2 | |
| (1) | (2) | |
| x1 | 1.98*** | 1.96*** |
| (0.03) | (0.03) | |
| Constant | 1.03*** | 10.97*** |
| (0.03) | (0.03) | |
| Observations | 1,000 | 1,000 |
| R2 | 0.80 | 0.79 |
| Adjusted R2 | 0.80 | 0.79 |
| Residual Std. Error (df = 998) | 1.00 | 1.03 |
| F Statistic (df = 1; 998) | 3,930.37*** | 3,649.38*** |
| Note: | p<0.1; p<0.05; p<0.01 | |
傾きはほとんど同じ: \(1.98\) と \(1.96\)
切片は異なる: \(1.03\) と \(10.97\)
二つのモデルを可視化してみる
シミュレーションの結果
・「誤差項の期待値がゼロでない」モデル \(y_2\) の場合でも、影響を受けるのは \(Y\) 切片 (1.03→10.97) だけであり、傾き
\(\beta_i = 1.98\)
にはほとんど影響はない
・誤差項 \(ε_i\)
の期待値はゼロと見なすことができる
回帰モデルにおけるパラメータ(母数)が線形という仮定
「線形回帰モデルでは、誤差を除く応答変数の決定要因が各説明変数の線形関数で表されなければならない」という仮定
この仮定は、最小二乗法による回帰係数の推定量が「不偏」であるために必要
応答変数を \(y\)、説明変数を \(x_m (m = 1, 2, ..., k )\)とすると、次の関係があるということ
\[y = β_1x_1 + β_2x_2 +...+ β_kx_k\]
\[y = x_1x_2x_3\]
★ 解決策 1 → 両辺の対数をとる
\[logy=log(x_1x_2x_3)=logx_1 +logx_2 +logx_3\]
\[Y = X_1 + X_2 + X_3\]
★ 解決策 2 → \(z = x_1x_2x_3\) と定義
\[y = z\]
\[y=logx\]
★ 解決策 → 新たに \(t=logx\) という変数を作る
log関数で変数変換し、\(y=t\) という関係を考えることができる| モデル | 解釈 |
|---|---|
| I. \(\hat{y}_i = \hat{\beta_0} + \hat{\beta_1}x_{1i}\) | \(x_1\)が1単位変化すると、\(y\) は\(\hat{\beta_1}単位変化する\) |
| II. \(\hat{y}_i = \hat{\beta_0} + \hat{\beta_1}log(x_{1i})\) | \(x_1\)が1%変化すると、\(y\) は\(\hat{\beta_1}/100単位変化する\) |
| III. \(\widehat{log(y_i)} = \hat{\beta_0} + \hat{\beta_1}x_{1i}\) | \(x_1\)が1単位変化すると、\(y\) は100\(\hat{\beta_1}\)%変化する |
| IV. \(\widehat{log(y_i)} = \hat{\beta_0} + \hat{\beta_1}log(x_{1i})\) | \(x_1\)が1%変化すると、\(y\) は\(\hat{\beta_1}\)%変化する |
gapminder
パッケージに含まれている「GDP」と「寿命」変数を使うgapminder)gapminder パッケージをダウンロードするlibrary(gapminder)gapminder
には次の変数が含まれているnames(gapminder)[1] "country" "continent" "year" "lifeExp" "pop" "gdpPercap"
lifeExp,
gdpPercapgapminder_1 <- gapminder |>
dplyr::select(year, lifeExp, gdpPercap) log_gdp,
log_lifeExp) を作り、変数を限定するgapminder_1 <- gapminder_1 |>
mutate(log_gdpPercap = log(gdpPercap)) |>
mutate(log_lifeExp = log(lifeExp)) |>
round(digits = 2) |> # 表示を 2 桁に指定
dplyr::filter(year == 2007) # 2007年のデータだけを残すgapminder の記述統計を示すlibrary(stargazer)| 変数名 | 詳細 |
|---|---|
| lifeExp | 寿命 |
| log_lifeExp | 寿命(対数変換) |
| gdpPercap | 1 人あたり GDP, 2005年の米ドル基準 |
| log_gdpPercap | 1 人あたり GDP(対数変化) |
{r, results = "asis"}
と指定するstargazer(as.data.frame(gapminder_1),
type ="html",
digits = 2)| Statistic | N | Mean | St. Dev. | Min | Max |
| year | 142 | 2,007.00 | 0.00 | 2,007 | 2,007 |
| lifeExp | 142 | 67.01 | 12.07 | 39.61 | 82.60 |
| gdpPercap | 142 | 11,680.07 | 12,859.94 | 277.55 | 49,357.19 |
| log_gdpPercap | 142 | 8.62 | 1.36 | 5.63 | 10.81 |
| log_lifeExp | 142 | 4.19 | 0.20 | 3.68 | 4.41 |
gapminder の中身を確認する
DT::datatable(gapminder_1)gdpPercap & lifeExpgdpPercap)」を x 軸、「平均寿命
(lifeExp)」を y 軸に指定して散布図を描いてみるgapminder_1 |>
ggplot() +
geom_point(aes(x = gdpPercap,
y = lifeExp)) +
labs(x = "1 人あたりGDP (USD)",
y = "寿命") +
theme_bw(base_family = "HiraKakuProN-W3")model_1 <- lm(lifeExp ~ gdpPercap, data = gapminder_1)model_4 と一緒に後に示すlog_gdpPercap & log_lifeExpGDP:
gdpPercap」の記述統計を見ると、データの範囲が
277米ドル(最低額)から 49,359
米ドル(最高額)まで非常に広範囲に及ぶ対数変換 した変数
(log_gdpPercap) と「平均寿命」(log_lifeExp)
との散布図を描いてみるgapminder_1 |>
ggplot() +
geom_point(aes(x = log_gdpPercap,
y = log_lifeExp)) +
labs(x = "1 人あたりGDP (USD)の対数値",
y = "寿命(対数変換)") +
theme_bw(base_family = "HiraKakuProN-W3")model_2 <- lm(log_lifeExp ~ log_gdpPercap, data = gapminder_1)model_1 と model_2
の分析結果を表示させるstargazer(model_1, model_2, type = "html")| Dependent variable: | ||
| lifeExp | log_lifeExp | |
| (1) | (2) | |
| gdpPercap | 0.001*** | |
| (0.0001) | ||
| log_gdpPercap | 0.113*** | |
| (0.008) | ||
| Constant | 59.566*** | 3.211*** |
| (1.010) | (0.067) | |
| Observations | 142 | 142 |
| R2 | 0.461 | 0.607 |
| Adjusted R2 | 0.457 | 0.605 |
| Residual Std. Error (df = 140) | 8.898 | 0.124 |
| F Statistic (df = 1; 140) | 119.556*** | 216.516*** |
| Note: | p<0.1; p<0.05; p<0.01 | |
対数変換した回帰分析 (model_4)
の解釈
・一人当たりGDPが 1ドル増えると、寿命が 0.001歳増える
→ 一人当たりGDPが 1000ドル増えると、寿命が 1 歳増える
\[Y_i = \beta_0 + \beta_1x_{1i} + .... + \beta_kx_{ki} + ε_i\]
「誤差が独立」というのは「別々の個体 \(i\) と \(j\) の誤差の間に何の関係もない」という意味
上に書いた回帰モデルでは、応答変数のうち説明変数に含まれていない変数の影響は全て誤差項に含まれる
従って、誤差が独立であるというのは「説明変数に含まれるもの以外で応答変数に影響する要因が互いに独立である(=無関係である)」ということ
「誤差が独立である」ことを式で表現すると
\[Cov (ε_i, ε_j) = 0\]
\[i ≠ j であるようなすべての (i , j) について、ε_i と ε_j の共分散が 0 になる\]
これは「系列相関がない」と同じ意味
例えば、個人のビールの消費量 \(Beer\) が個人の収入 \(Income\) によって影響を受けるという次のようなモデルを考える
\[Beer_i = β_0 - β_1Income_i + ε_i\]
「誤差が独立である」=「2
つの個体、\(i\) と
\(j\)
の誤差の間に何の関係もない」
上記の回帰モデルでは、応答変数 \(Beer\) のうち説明変数 \(Income\) に含まれていない変数の影響は誤差項 \(ε_i\)に含まれる
従って、「誤差が独立である」=「説明変数に含まれるもの以外で、応答変数に影響する要因が互いに独立である(=無関係である)」
☆ \(i\) さんと \(j\)
さんはゼミの飲み友達で、一緒にビールを飲みに行く
→ \(i\)
さんのビール消費量が増えると、\(j\)
さんのビール消費量も増える
「平均独立」と「条件付き平均独立」
mean independent) であるというconditional mean independent) であるという「誤差が独立」であるという仮定は最小二乗法による回帰係数の推定量が不偏であるために必要
→ 「統制すべき交絡因子が十分にモデルに含まれていなければならない」ということにつながる
どのような変数を交絡因子としてモデルに含めるべきか(含めるべきではないか)という問題は次のセッションを参照:
「26.
重回帰分析 9(コントロール変数の選び方:バックドア基準と
DAG)」
不要な変数をモデルに含めた場合と、必要な変数をモデルに含めなかった場合にどういうことが起きるかという問題は次のセッションを参照:
「25.
重回帰分析(脱落変数バイアス・処置後変数バイアス)」
\[全ての個体 i について、Var(ε_i) = σ^2\]
homoskerdasticity |
均一分散 | 分散が均一であること |
heteroskedasticity |
不均一分散 | 分散が均一でないこと |
均一分散(左側の図)と不均一性分散(右側の図)の事例
heteroskerdasticity)
の事例としては「月収」と「食費」が考えられる社会科学の事例における説明変数と誤差の分散
応答変数のばらつきが大きくなるほど
→ 誤差のばらつきが大きくなる
→ 説明変数で説明できない部分が大きくなる
→ 誤差に分散不均一性が生じる
→ 標準誤差に影響がでて、統計的な有意性が得にくくなる
・ 誤差項の分散が不均一でも、最小二乗法による回帰係数の推定量に偏りは発生しない
Breusch-Pagan Test
によって、誤差項の分散が均一かどうか検定できる
政府のパフォーマンスに関するデータ
(putnam.csv) をダウンロードし、RProject フォルダー内の
data フォルダに保存する
データのダウンロードが終わったら、データを読み込み
putnam と名前を付ける
putnam <- read_csv("data/putnam.csv")gov_p を応答変数とした単回帰分析を行うfit.putnam <- lm(gov_p ~ cc, data = putnam) residuals) をプロットするres.putnam <- resid(fit.putnam) # 残差を取り出す
pred.putnam <- fitted(fit.putnam) # 予測値を取り出す
plot(pred.putnam,
res.putnam,
main = "Residuals vs Fitted",
xlab = "Fitted Values",
ylab = "residuals")
abline(h = 0, col = "red") # 残差 = 0 に平行な赤線を引く応答変数 Fitted Values
の大小によって、誤差項が偏っている(=
不均一分散)ようには見えない
Rパッケージ lmtest の
bytest()関数を使った不均一分散の診断
library(lmtest)bptest(fit.putnam)
studentized Breusch-Pagan test
data: fit.putnam
BP = 0.36787, df = 1, p-value = 0.5442
bytest()関数の診断基準p-valueが 0.05以上の場合 |
誤差項の分散が均一である | |
p-valueが 0.05以下の場合 |
誤差項の分散が不均一である | |
p-value = 0.5442 で 0.05 以上なので、「誤差項の分散が均一である」Breusch-Pagan Test
によって、誤差項の分散が均一かどうか検定できる
データのダウンロード方法
RProject フォルダ内に data
という名称のフォルダを作成する
hr96-21.csv をクリックしてデータをパソコンにダウンロード
RProject フォルダ内に data
という名称のフォルダを作成する
ダウンロードした hr96-21.csv
を手動でRProject フォルダ内にある data
フォルダに入れる
選挙データの読み取り方法
次のいずれかの方法で
hr96-21.csv を読み取る
na = "."というコマンドは「欠損値をドットで置き換える」という意味numeric)」型のデータが「」文字型
(character)」として認識されるなど、エラーの原因になるため、読み取る時点で事前に対処するhr <- read_csv("data/hr96-21.csv",
na = ".") locale()関数を使って日本語エンコーディング
(cp932) を指定するhr <- read.csv("data/hr96-21.csv",
na = ".",
locale = locale(encoding = "cp932"))hr <- read.csv("data/hr96-21.csv",
na = ".") names(hr) [1] "year" "pref" "ku" "kun"
[5] "wl" "rank" "nocand" "seito"
[9] "j_name" "gender" "name" "previous"
[13] "age" "exp" "status" "vote"
[17] "voteshare" "eligible" "turnout" "seshu_dummy"
[21] "jiban_seshu" "nojiban_seshu"
hr09 <- hr |>
dplyr::filter(year == 2009) |>
dplyr::select(voteshare, age)voteshare を応答変数とした単回帰分析を行うfit.hr09 <- lm(voteshare ~ age, data = hr09) residuals)
のプロットを描いてみるfit.hr09 <- lm(voteshare ~ age, data = hr09)
res.hr09 <- resid(fit.hr09) # 残差を取り出す
pred.hr09 <- fitted(fit.hr09) # 予測値を取り出す
plot(pred.hr09,
res.hr09,
main = "Residuals vs Fitted",
xlab = "Fitted Values",
ylab = "residuals")
abline(h = 0, col = "red") # 残差 = 0 に平行な赤線を引くFitted Values
が大きくなるにつれて、Risiduals(赤の平行点線)を中心に下の方がデータが広がっているので、誤差の独立性は保たれていないようにみえる
Rパッケージ lmtest の
bytest()関数を使った不均一分散の診断
library(lmtest)bptest(fit.hr09)
studentized Breusch-Pagan test
data: fit.hr09
BP = 13.982, df = 1, p-value = 0.0001846
bytest()関数の診断基準p-valueが 0.05以上の場合 |
誤差項の分散が均一である | |
p-valueが 0.05以下の場合 |
誤差項の分散が不均一である | |
p-value = 0.00509 で 0.05 以下なので、「誤差項の分散が不均一である」| 1. 加重最小二乗法 | 重み付けたOLS | 不均一にしている変数とその関数形が分かる場合 |
| 2. 不均一分散に頑健な標準誤差 | 標準誤差を修正 | 不均一にしている変数とその関数形が分からない場合 |
WLS: Weighted least squares)precisionと呼ぶ)に比例し
た重みを使った「重み付き最小二乗法」の推定量を使う2. 誤差項の期待値がゼロ
6. 誤差の分散が均一である
\[ε_i 〜 N (0, σ^2)\]
残差が正規分布に従っていないケース(左)と正規分布に従っているケース(右)の事例
BLUE: Best Linear Unbiased Estimator)である 5
つの仮定の中に「誤差が正規分布に従う」という仮定は含まれていないBLUE) であるJarque-Bera test for normality)normtest の
jp.norm.test()関数を使って検定可能bytest()関数の診断基準p-valueが 0.05以上の場合 |
誤差項が正規分布に従っている |
p-valueが 0.05以下の場合 |
誤差項が正規分布に従っていない |
putnam <- read_csv("data/putnam.csv")
fit.putnam <- lm(gov_p ~ cc, data = putnam)
res.putnam <- resid(fit.putnam) # 残差を取り出す
par(family = "HiraKakuProN-W3") # 日本語文字化け防止 (Macユーザのみ)
qqnorm(res.putnam, main = "残差の正規QQプロット",
xlab = "標準正規分布", ylab = "標準化した残差の分布")
qqline(res.putnam, col = "red") # 正規分布に一致する場合の直線を引くデータが直線の周りに偏りなく分布しているので、正規QQプロットによる誤差が正規分布に従っているといえる
2009総選挙の残差のヒストグラムを描いてみる
hist(res.putnam)bytest()関数を使って正規性を診断をするlibrary(normtest)resid_putnam <- residuals(fit.putnam)
jb.norm.test(resid_putnam)
Jarque-Bera test for normality
data: resid_putnam
JB = 1.3857, p-value = 0.2475
p-value = 0.2475
なので帰無仮説「誤差項が正規分布に従っている」は棄却できないhr09 <- read_csv("data/hr96-21.csv") |>
filter(year == 2009)
fit.hr09 <- lm(voteshare ~ age, data = hr09)
res.hr09 <- resid(fit.hr09) # 残差を取り出す
par(family = "HiraKakuProN-W3") # 日本語文字化け防止 (Macユーザのみ)
qqnorm(res.hr09, main = "残差の正規QQプロット",
xlab = "標準正規分布", ylab = "標準化した残差の分布")
qqline(res.hr09, col = "red") # 正規分布に一致する場合の直線を引くデータが直線の周りに偏りなく分布しているとはいえない
2009総選挙の残差のヒストグラムを描いてみる
hist(res.hr09)bytest()関数を使って正規性を診断をするlibrary(normtest)resid_hr09 <- residuals(fit.hr09)
jb.norm.test(res.hr09)
Jarque-Bera test for normality
data: res.hr09
JB = 73.299, p-value < 2.2e-16
p-value < 2.2e-16
なので帰無仮説「誤差項が正規分布に従っている」は棄却par( ) 関数を使って 2 行 2
列のグラフィック環境を整えてから、plot( )関数を使って、4
つのプロットを並べて表示することもできるfit.putnam <- lm(gov_p ~ cc + econ, data = putnam) # gov_p を応答変数とした重回帰分析を行う par(mfrow = c(2,2)) # 2 行 2 列のグラフィック環境を作る plot(fit.putnam) # 4 つのプロットを表示
多重共線性 (multicollinearity) の定義:
「説明変数の間の相関係数が 1
ではないものの、強い相関があること」
次のようなモデルを考える
\[Y_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + ε_i\]
出典:高橋(2022)p.107を参考に著者が作成した図
VIFVIF: variance inflation factor)
という指標を使って診断する\[\hat{Y} = \hat{\beta_0} + \hat{\beta_1X_{1i}} + \hat{\beta_2X_{2i}} +...+ \hat{\beta_pX_{pi}}\]
\[\hat{X_{ji}} = \hat{γ_0} + \hatγX_{1i} ... + \hatγX_ {pi} \]
\[VIF = \frac{1}{1-R^2_j}\]
VIFRパッケージ carを使ってモデルの VIF
を計算できるputnam <- read_csv("data/putnam.csv")
fit.putnam <- lm(gov_p ~ cc + econ,
data = putnam) # gov_p を応答変数とした重回帰分析を行うlibrary("car")
vif(fit.putnam) cc econ
5.632172 5.632172
VIF > 10 なら、多重共線性の疑いがあるVIF > 10 というのは説明変数間の相関係数が
0.9 を超えるような場合cc と econ の VIF は 5.6
程度なので問題ないVIF > 10 の時の対処法VIFage
が小選挙区でのランク rank に与える影響rankageage
が最も興味のある重要な変数| 変数名 | 詳細 | |
|---|---|---|
rank |
小選挙区でのランク | 応答変数 |
age |
候補者の年齢 | 説明変数 |
vote |
候補者が得た票数 | 共変量 |
voteshare |
候補者が得た得票率 (%) | 共変量 |
vote」と「候補者が得た得票率
voteshare」をモデルに含めるvote も voteshare
も「交絡因子」ではないrank」に影響を与えるし、年齢とは「相関」があると思われるため便宜上モデルに含めるhr <- read_csv("data/hr96-21.csv")hr09 <- hr |>
dplyr::filter(year == 2009) |> # 2009年総選挙に限る
dplyr::select(age, vote, voteshare, rank) # 必要な変数を選ぶstargazer(as.data.frame(hr09),
type = "html",
digits = 2)| Statistic | N | Mean | St. Dev. | Min | Max |
| vote | 1,139 | 61,939.59 | 53,880.34 | 177 | 201,461 |
| voteshare | 1,139 | 26.34 | 22.20 | 0.10 | 95.30 |
| rank | 1,139 | 2.50 | 1.24 | 1 | 9 |
hr09 |>
ggplot(aes(age, as.factor(rank))) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw()hr09 |>
ggplot(aes(age, vote)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw()hr09$age <- as.numeric(hr09$age)
hr09$vote <- as.numeric(hr09$vote)with(hr09, cor.test(age, vote))
Pearson's product-moment correlation
data: age and vote
t = 6.5822, df = 1133, p-value = 7.072e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1352367 0.2473405
sample estimates:
cor
0.1919146
hr09 |>
ggplot(aes(vote, voteshare)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw()with(hr09, cor.test(vote, voteshare))
Pearson's product-moment correlation
data: vote and voteshare
t = 117.15, df = 1137, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.9562771 0.9651931
sample estimates:
cor
0.960984
m1 <- lm(rank ~ age, data = hr09)
m2 <- lm(rank ~ age + vote, data = hr09)
m3 <- lm(rank ~ age + voteshare, data = hr09)
m4 <- lm(rank ~ age + vote + voteshare, data = hr09)stargazer(m1, m2, m3, m4,
type = "html")| Dependent variable: | ||||
| rank | ||||
| (1) | (2) | (3) | (4) | |
| age | -0.021*** | -0.003 | -0.001 | -0.001 |
| (0.003) | (0.002) | (0.002) | (0.002) | |
| vote | -0.00002*** | -0.00000 | ||
| (0.00000) | (0.00000) | |||
| voteshare | -0.050*** | -0.046*** | ||
| (0.001) | (0.003) | |||
| Constant | 3.538*** | 3.847*** | 3.845*** | 3.847*** |
| (0.166) | (0.086) | (0.077) | (0.077) | |
| Observations | 1,135 | 1,135 | 1,135 | 1,135 |
| R2 | 0.035 | 0.743 | 0.793 | 0.794 |
| Adjusted R2 | 0.034 | 0.743 | 0.793 | 0.793 |
| Residual Std. Error | 1.221 (df = 1133) | 0.630 (df = 1132) | 0.566 (df = 1132) | 0.566 (df = 1131) |
| F Statistic | 41.288*** (df = 1; 1133) | 1,637.589*** (df = 2; 1132) | 2,168.995*** (df = 2; 1132) | 1,448.713*** (df = 3; 1131) |
| Note: | p<0.1; p<0.05; p<0.01 | |||
age」が候補者の「ランク」に与える影響vote」と「候補者が得た得票率
voteshare」の偏回帰係数は気にする必要なしage」偏回帰係数Model_1:age が 1 歳増えると、候補者のランクが 0.02
ポイントだけ良くなる(= 1 位に近くなる)Model_1
には統制変数が含まれていないので、この結果は候補者の年齢
age を過大評価していると思われるModel_2 &
Model_3) を推定する必要ありModel_2 & Model_3vote をモデルに含めると、候補者の年齢
age の影響力は消えるvoteshare
をモデルに含めると、候補者の年齢 age の影響力は消えるModel_1 の結果は候補者の年齢 age
を課題評価していることが裏付けられたModel_4vote
と「候補者の得票率」voteshare
を両方同時にモデルに含めても、候補者の年齢 age
の影響力は消え、しかも、Model_2, 3, 4と、候補者の年齢
ageの係数が変わらないage
の影響力はない」という推定が正しそうModel_4 の VIF を計算するcarパッケージを使って、モデル4の VIF
を計算するlibrary(car)
vif(m4) age vote voteshare
1.043137 13.127862 13.187917
多重共線のポイント
・ 回帰モデルが含む 2 つの共変量 vote と
voteshare の VIF が 13 以上
・共変量の間に多重共線が起きている
→ VIF が 10 以上ではあるが、これら
2 つの変数をモデルに加えても「候補者の年齢
age」は正しく推定できているので問題なし
・共変量の間に多重共線が起きていても、特に問題視する必要はない
・むしろ必要な共変量をモデルに含めないと、交絡を取り除けない
・必要と思われる共変量はモデルに含めるべき
heteroskedasticity robust standard error)を考慮した回帰分析を行うheteroskedasticity robust standard error)を考慮した回帰分析を行うlm()
の中で、個人や時点などをダミー変数としてモデルに含めるgapminder
パッケージに含まれている「GDP」と「寿命」変数を使うデータの準備 (gapminder)
gapminder パッケージをダウンロードするlibrary(gapminder)gapminder
には次の変数が含まれているnames(gapminder)[1] "country" "continent" "year" "lifeExp" "pop" "gdpPercap"
lifeExp,
gdpPercapgapminder <- gapminder |>
dplyr::select(year, lifeExp, gdpPercap) log_gdp,
log_lifeExp) を作り、変数を限定するgapminder_2 <- gapminder |>
mutate(log_gdpPercap = log(gdpPercap)) |>
mutate(log_lifeExp = log(lifeExp)) |>
round(digits = 2) # 表示を 2 桁に指定gapminder の記述統計を示すlibrary(stargazer)| 変数名 | 詳細 |
|---|---|
| year | 年度 |
| lifeExp | 寿命 |
| log_lifeExp | 寿命(対数変換) |
| gdpPercap | 1 人あたり GDP, 2005年の米ドル基準 |
| log_gdpPercap | 1 人あたり GDP(対数変化) |
{r, results = "asis"}
と指定するstargazer(as.data.frame(gapminder_2),
type ="html",
digits = 2)| Statistic | N | Mean | St. Dev. | Min | Max |
| year | 1,704 | 1,979.50 | 17.27 | 1,952 | 2,007 |
| lifeExp | 1,704 | 59.47 | 12.92 | 23.60 | 82.60 |
| gdpPercap | 1,704 | 7,215.33 | 9,857.45 | 241.17 | 113,523.10 |
| log_gdpPercap | 1,704 | 8.16 | 1.24 | 5.49 | 11.64 |
| log_lifeExp | 1,704 | 4.06 | 0.23 | 3.16 | 4.41 |
gapminder の中身を確認する
DT::datatable(gapminder_2)gdpPercap & lifeExpgdpPercap)」を x 軸、「平均寿命
(lifeExp)」を y 軸に指定して散布図を描いてみるgapminder_2 |>
ggplot() +
geom_point(aes(x = gdpPercap,
y = lifeExp)) +
labs(x = "1 人あたりGDP (USD)",
y = "寿命") +
theme_bw(base_family = "HiraKakuProN-W3")log_gdpPercap & log_lifeExpGDP:
gdpPercap」の記述統計を見ると、データの範囲が
241米ドル(最低額)から 113,523
米ドル(最高額)まで非常に広範囲に及ぶ対数変換 した変数
(log_gdpPercap) と「平均寿命」(log_lifeExp)
との散布図を描いてみるgapminder_2 |>
ggplot() +
geom_point(aes(x = log_gdpPercap,
y = log_lifeExp)) +
labs(x = "1 人あたりGDP (USD)の対数値",
y = "寿命(対数変換)") +
theme_bw(base_family = "HiraKakuProN-W3")model_22 <- lm(log_lifeExp ~ log_gdpPercap, data = gapminder_2)model_22 の分析結果を表示させるstargazer(model_22, type = "html")| Dependent variable: | |
| log_lifeExp | |
| log_gdpPercap | 0.147*** |
| (0.003) | |
| Constant | 2.864*** |
| (0.023) | |
| Observations | 1,704 |
| R2 | 0.613 |
| Adjusted R2 | 0.613 |
| Residual Std. Error | 0.145 (df = 1702) |
| F Statistic | 2,698.395*** (df = 1; 1702) |
| Note: | p<0.1; p<0.05; p<0.01 |
対数変換した回帰分析
(model_22) の解釈
year ダミー を含めるfactor(year) をモデルに追加するmodel_33 <- lm(log_lifeExp ~ log_gdpPercap + factor(year),
data = gapminder_2)stargazer(model_22, model_33,
type = "html")| Dependent variable: | ||
| log_lifeExp | ||
| (1) | (2) | |
| log_gdpPercap | 0.147*** | 0.135*** |
| (0.003) | (0.003) | |
| factor(year)1957 | 0.033** | |
| (0.016) | ||
| factor(year)1962 | 0.058*** | |
| (0.016) | ||
| factor(year)1967 | 0.079*** | |
| (0.016) | ||
| factor(year)1972 | 0.095*** | |
| (0.016) | ||
| factor(year)1977 | 0.116*** | |
| (0.016) | ||
| factor(year)1982 | 0.144*** | |
| (0.016) | ||
| factor(year)1987 | 0.171*** | |
| (0.016) | ||
| factor(year)1992 | 0.184*** | |
| (0.016) | ||
| factor(year)1997 | 0.186*** | |
| (0.016) | ||
| factor(year)2002 | 0.185*** | |
| (0.016) | ||
| factor(year)2007 | 0.185*** | |
| (0.016) | ||
| Constant | 2.864*** | 2.841*** |
| (0.023) | (0.023) | |
| Observations | 1,704 | 1,704 |
| R2 | 0.613 | 0.684 |
| Adjusted R2 | 0.613 | 0.681 |
| Residual Std. Error | 0.145 (df = 1702) | 0.131 (df = 1691) |
| F Statistic | 2,698.395*** (df = 1; 1702) | 304.605*** (df = 12; 1691) |
| Note: | p<0.1; p<0.05; p<0.01 | |
year の値を確認してみるtable(gapminder_2$year)
1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 2002 2007
142 142 142 142 142 142 142 142 142 142 142 142
factor(year)1957 から
factor(year)2007 まで 11 個しかない0.0330.033 だけ増える」と解釈する固定効果を含むモデル(model_33)
の解釈
estimatr パッケージの
lm_robust()関数を使っても可能library(estimatr)fixed_robust <- lm_robust(log_lifeExp ~ log_gdpPercap,
data = gapminder_2,
fixed_effects = ~ year,
se_type = "classical")summary(fixed_robust)
Call:
lm_robust(formula = log_lifeExp ~ log_gdpPercap, data = gapminder_2,
fixed_effects = ~year, se_type = "classical")
Standard error type: classical
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
log_gdpPercap 0.1347 0.002639 51.06 0 0.1295 0.1399 1691
Multiple R-squared: 0.6837 , Adjusted R-squared: 0.6815
Multiple R-squared (proj. model): 0.6065 , Adjusted R-squared (proj. model): 0.6037
F-statistic (proj. model): 2607 on 1 and 1691 DF, p-value: < 2.2e-16